In this problem, we are looking for a scalar pressure field
\( {\cal R}^3 \ni (x,y,z)\times [0,T] \rightarrow u(x,y,z;t) \in {\cal R } \)
in a heterogeneous area. This problem is related to the issue of oil extraction. We assume that the area in which we model our flow consists of heterogeneous material representing oil deposits located in the crevices of geological layers. We assume that boreholes were made in which a pump pumping liquid under high pressure was placed, the purpose of which is to pump out the oil reservoir, in particular towards suction pumps located in separate wells.
Equation modeling pressure changes due to pump operation
\( \frac{\partial u}{\partial t } - \nabla \cdot (\kappa \exp(\mu u)\,\nabla u) = f \)
which after rewritting
\( \frac{\partial u(x,y,z;t)}{\partial t } - \frac{\partial }{\partial x} \left(\kappa(x,y,z) \exp(\mu(x,y,z) u(x,y,z;t)) \frac{\partial u(x,y,z;t)}{\partial x} \right) + \\ + \frac{\partial }{\partial y} \left(\kappa(x,y,z; u(x,y,z;t)) \exp(\mu u(x,y,z;t)) \frac{\partial u(x,y,z;t)}{\partial y } \right) = f(x,y,z;t) \)
In the equation above, the function \( \kappa(x,y,z) \) determines the permeability of the area. This function is shown in Fig. 1. The value range varies from 0 to 1000. Value 0 means no permeability, value 1000 means very high permeability (e.g. a crude oil fracture in a rock).

Factor \( \mu=1000 \) is a geological model parameter.
Note that the nonlinear flow equation is very similar to the heat transport equation. The only difference is the non-linear coefficient
\( \kappa(x,y,z) \exp(\mu(x,y,z) u(x,y,z;t) \)
which means that the rate of pressure changes in a given area changes proportionally to the permeability coefficient and to the current pressure value, and it happens in a non-linear manner, determined by the exponential function. This coefficient is the cause of the non-linearity of our equation. If this term was left on the left-hand side of the equation, then this equation would require linearization, for example, using the Newton-Raphson algorithm. In the case of nonlinear problems, however, it is possible to use the explicit method which gives the possibility to shift the non-linear term to the right-hand side of the equation.
We approximate the derivative over time using the finite difference method
\( u_{t+1}-u_t - dt* \nabla \cdot (\kappa \exp(\mu u)\,\nabla u) = dt*f \)
The right-hand side function models the pumping and suction pumps
\( f(x, t) = \textcolor{blue}{\sum_{p \in P} \phi\left(\left\|x_p - x\right\|\right)} - \textcolor{red}{\sum_{s \in S} u(x, t) \phi\left(\left\|x_s - x\right\|\right) } \)
Where \( \textcolor{blue}{\sum_{p \in P} \phi\left(\left\|x_p - x\right\|\right) } \) are pumping pumps (providing non-zero pressure) and \( -\textcolor{red}{\sum_{s \in S} \phi\left(\left\|x_p - x\right\|\right) } \) are suction pumps ("delivering" negative pressure). The sum goes over all positions of the P pumping and S suction pumps in our modeled system.
The pressure applied by the pump drops with the distance from the center of the pump according to the formula also illustrated in Fig. 2
\( \phi(t) = \begin{cases} \left(\frac{t}{r} - 1\right)^2 \left(\frac{t}{r} + 1\right)^2 & \text{for } t \leq r \\ 0 & \text{for } t > r \\ \end{cases} \).

In order to simulate using the isogeometric finite element method, we introduce discretization using the B-spline basis functions
\( \{B^x_{i,2}(x)B^y_{j,2}(y)B^z_{k,2}(z)\}_{i=1,...,N_x;j=1,...,N_y;k=1,...,N_z } \)
They will be used to approximate the simulated scalar pressure field to the current time instant
\( u(x,y,z;t+1)\approx \sum_{i=1,...,N_x;j=1,...,N_y;k=1,...,N_z} u^{t+1}_{i,j,k } B^x_i(x)B^y_j(y)B^z_k(z) \)
Similarly, for testing, we will use the B-spline basis functions:
\( \{ B^x_l(x)B^y_m(y)B^z_n(z) \}_{l=1,...,N_x;m=1,...,N_y;n=1,...,N_z } \)
Our equation is multiplied by the test function and integrated over the area of the three-dimensional simulation (we mark the integration with a parenthesis representing the dot product in the L2 space)
\( (u_{t+1},w)=(u_t,w) - dt ( \nabla \cdot (\kappa \exp(\mu u)\,\nabla u),w)+dt (f,w) \)
We can isolate the term on the right-hand side through the parts, assuming that our area is so large that the pressure changes on the edge are zero (then the marginal integral resulting from integration by parts will disappear)
\( (u_{t+1},w)=(u_t,w) + dt ( \kappa \exp(\mu u)\,\nabla u,\nabla w)+dt (f,w) \).

As in the general case of the explicit method, on the left-hand side we now have a mass matrix that can be factored in a linear time, and on the right-hand side we have our "bitmap" which describes the state at the previous moment plus changes during the time step dt resulting from the physics of nonlinear flow in a heterogeneous medium and modification of pressure fields in the vicinity of pumping and suction pumps.